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We here discuss the emergence of Quasi Stationary States (QSS), a universal feature of systems 
with long-range interactions. With reference to the Hamiltonian Mean Field (HMF) model, numeri- 
cal simulations are performed based on both the original iV-body setting and the continuum Vlasov 
model which is supposed to hold in the thermodynamic limit. A detailed comparison unambiguously 
demonstrates that the Vlasov-wave system provides the correct framework to address the study of 
QSS. Further, analytical calculations based on Lynden-Bell's theory of violent relaxation are shown 
to result in accurate predictions. Finally, in specific regions of parameters space, Vlasov numerical 
solutions are shown to be affected by small scale fluctuations, a finding that points to the need for 
novel schemes able to account for particles correlations. 
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The Vlasov equation constitutes a universal theoreti- 
cal framework and plays a role of paramount importance 
in many branches of applied and fundamental physics. 
Structure formation in the universe is for instance a 
rich and fascinating problem of classical physics: The 
fossile radiation that permeates the cosmos is a relic 
of microfluctuation in the matter created by the Big 
Bang, and such a small perturbation is believed to have 
evolved via gravitational instability to the pronounced 
agglomerations that we see nowdays on the galaxy clus- 
ter scale. Within this scenario, gravity is hence the en- 
gine of growth and the Vlasov equation governs the dy- 
namics of the non baryonic "dark matter" Further- 
more, the continuous Vlasov description is the reference 
model for several space and laboratory plasma applica- 
tions, including many interesting regimes, among which 
the interpretation of coherent electrostatic structures ob- 
served in plasmas far from thermodynamic equilibrium. 
The Vlasov equation is obtained as the mean-field limit 
of the iV-body Liouville equation, assuming that each 
particle interacts with an average field generated by all 
plasma particles (i.e. the mean electromagnetic field de- 
termined by the Poisson or Maxwell equations where the 
charge and current densities are calculated from the par- 
ticle distribution function) while inter-particle correla- 
tions are completely neglected. 

Numerical simulations are presently one of the most 
powerful resource to address the study of the Vlasov 
equation. In the plasma context, the Lagrangian 
Part icle-In- Cell approach is by far the most popular, 
while Eulerian Vlasov codes are particularly suited for 
analyzing specific model problems, due to the associated 
low noise level which is secured even in the non-linear 



regime [2[. However, any numerical scheme designed to 
integrate the continuous Vlasov system involves a dis- 
cretization over a finite mesh. This is indeed an unavoid- 
able step which in turn affects numerical accuracy. A nu- 
merical (diffusive and dispersive) characteristic length is 
in fact introduced being at best comparable with the grid 
mesh size: as soon as the latter matches the typical length 
scale relative to the (dynamically generated) fluctuations 
a violation of the continuous Hamiltonian character of the 
equations occurs (see Refs. [3[). It is important to em- 
phasize that even if such non Vlasov effects are strongly 
localized (in phase space), the induced large scale topo- 
logical changes will eventually affect the system globally. 
Therefore, aiming at clarifying the problem of the valid- 
ity of Vlasov numerical models, it is crucial to compare a 
continuous Vlasov, but numerically discretized, approach 
to a homologous N-body model. 

Vlasov equation has been also invoked as a reference 
model in many interesting one dimensional problems, and 
recurrently applied to the study of wave-particles inter- 
acting systems. The Hamiltonian Mean Field (HMF) 
model [4], describing the coupled motion of TV rotators, 
is in particular assimilated to a Vlasov dynamics in the 
thermodynamic limit on the basis of rigorous results 
The HMF model has been historically introduced as rep- 
resenting gravitational and charged sheet models and is 
quite extensively analyzed as a paradigmatic representa- 
tive of the broader class of systems with long-range inter- 
actions [6[. A peculiar feature of the HMF model, shared 
also by other long-range interacting systems, is the pres- 
ence of Quasi Stationary States (QSS). During time evo- 
lution, the system gets trapped in such states, which are 
characterized by non Gaussian velocity distributions, be- 
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fore relaxing to the final Boltzmann-Gibbs equilibrium 
0] . An attempt has been made [1] to interpret the emer- 
gence of QSSs by invoking Tsallis statistics [9]. This 
approach has been later on criticized in [10], where QSSs 
were shown to correspond to stationary stable solutions 
of the Vlasov equation, for a particular choice of the ini- 
tial condition. More recently, an approximate analytical 
theory, based on the Vlasov equation, which derives the 
QSSs of the HMF model using a maximum entropy prin- 
ciple, was developed in [ll|. This theory is inspired by 
the pioneering work of Lynden-Bell [12[ and relies on 
previous work on 2D turbulence by Chavanis [Hj]. How- 
ever, the underlying Vlasov ansatz has not been directly 
examined and it is recently being debated [141 ]. 

In this Letter, we shall discuss numerical simulations 
of the continuous Vlasov model, the kinetic counterpart 
of the discrete HMF model. By comparing these results 
to both direct N-body simulations and analytical predic- 
tions, we shall reach the following conclusions: (i) the 
Vlasov formulation is indeed ruling the dynamics of the 
QSS; (ii) the proposed analytical treatment of the Vlasov 
equation is surprisingly accurate, despite the approxima- 
tions involved in the derivation; (iii) Vlasov simulations 
are to be handled with extreme caution when exploring 
specific regions of the parameters space. 

The HMF model is characterized by the following 
Hamiltonian 
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where 9j represents the orientation of the j-th rotor and 
Pj is its conjugate momentum. To monitor the evolution 
of the system, it is customary to introduce the magneti- 
zation, a macroscopic order parameter defined as M = 
|M| = |^m^|/iV, where = (cos 9^ sin 9i) stands 
for the microscopic magnetization vector. As previously 
reported [4|, after an initial transient, the system gets 
trapped into Quasi-Stationary States (QSSs), i.e. non- 
equilibrium dynamical regimes whose lifetime diverges 
when increasing the number of particles N . Importantly, 
when performing the mean-field limit (TV — > oo) before 
the infinite time limit, the system cannot relax towards 
Boltzmann-Gibbs equilibrium and remains permanently 
confined in the intermediate QSSs. As mentioned above, 
this phenomenology is widely observed for systems with 
long-range interactions, including galaxy dynamics [TBI, 
free electron lasers [TBI , 2D electron plasmas [TtJ . 

In the N — > oo limit the discrete HMF dynamics re- 
duces to the Vlasov equation 

df/dt + pdf/d6 -(a!V/d6)df/dp = 0, (2) 

where f(6,p,t) is the microscopic one-particle distribu- 
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The specific energy h[f] = J J (p 2 /2)f(0,p,t)d0dp - 
{Ml + M 2 - l)/2 and momentum P[f] 
j f pf(9,p,t)d9dp functionals are conserved quanti- 
ties. Homogeneous states are characterized by M = 0, 
while non- homogeneous states correspond to M ^ 0. 

Rigorous mathematical results [B[ demonstrate that, 
indeed, the Vlasov framework applies in the continuum 
description of mean-field type models. This observation 
corroborates the claim that any theoretical attempt to 
characterize the QSSs should resort to the above Vlasov 
based interpretative picture. Despite this, the QSS non- 
Gaussian velocity distributions have been fitted [8[ us- 
ing Tsallis' g-exponentials, and the Vlasov formalism as- 
sumed valid only for the limiting case of homogeneous 
initial conditions [T3 |. In a recent paper [TTj], the afore- 
mentioned velocity distribution functions were instead 
reproduced with an analytical expression derived from 
the Vlasov scenario, with no adjustable parameters and 
for a large class of initial conditions, including inhomo- 
geneous ones. The key idea dates back to the seminal 
work by Lynden-Bell [12] (see also 0, 0) and con- 
sists in coarse-graining the microscopic one-particle dis- 
tribution function f(9,p,t) by introducing a local aver- 
age in phase space. It is then possible to associate an 
entropy to the coarse-grained distribution /: The corre- 
sponding statistical equilibrium is hence determined by 
maximizing such an entropy, while imposing the conser- 
vation of the Vlasov dynamical invariants, namely en- 
ergy, momentum and norm of the distribution. We shall 
here limit our discussion to the case of an initial single 
particle distribution which takes only two distinct val- 
ues: /o = 1/(4A6>A P ), if the angles (velocities) lie within 
an interval centered around zero and of half-width A# 
(A p ), and zero otherwise. This choice corresponds to the 
so-called "water-bag" distribution which is fully specified 
by energy h[f] = e, momentum P[f] = o~ and the initial 
magnetization Mo = (M x o,M y o). The maximum en- 
tropy calculation is then performed analytically [TT[ and 
results in the following form of the QSS distribution 
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where /3//o, A//o and /i//o are rescaled Lagrange mul- 
tipliers, respectively associated to the energy, momen- 
tum and normalization. Inserting expression ([6]) into the 
above constraints and recalling the definition of M x [f], 
M y [f], one obtains an implicit system which can be 
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solved numerically to determine the Lagrange multipliers 
and the expected magnetization in the QSS. Note that 
the distribution ([6]) differs from the usual Boltzmann- 
Gibbs expression because of the "fermionic" denomina- 
tor. Numerically computed velocity distributions have 
been compared in [11] to the above theoretical predictions 
(where no free parameter is used), obtaining an overall 
good agreement. However, the central part of the distri- 
butions is modulated by the presence of two symmetric 
bumps, which are the signature of a collective dynamical 
phenomenon [11]. The presence of these bumps is not 
explained by our theory. Such discrepancies has been re- 
cently claimed to be an indirect proof of the fact that 
the Vlasov model holds only approximately true. We 
shall here demonstrate that this claim is not correct and 
that the deviations between theory and numerical obser- 
vation are uniquely due to the approximations built in 
the Lynden-Bell approach. 

A detailed analysis of the Lynden-Bell equilibrium 
(j6]) in the parameter plane (Mo,e) enabled us to un- 
ravel a rich phenomenology, including out of equilibrium 
phase transitions between homogeneous (Mqss — 0) 
and non-homogeneous (Mqss 7^ 0) QSS states. Sec- 
ond and first order transition lines are found that sep- 
arate homogeneous and non homogeneous states and 
merge into a tricritical point approximately located in 
(Mo,e) = (0.2,0.61). When the transition is second or- 
der two extrema of the Lynden-Bell entropy are identified 
in the inhomogeneous phase: the solution M = corre- 
sponds to a saddle point, being therefore unstable; the 
global maximum is instead associated to M / 0, which 
represents the equilibrium predicted by the theory. This 
argument is important for what will be discussed in the 
following. 

Let us now turn to direct simulations, with the aim of 
testing the above scenario, and focus first on the kinetic 
model ([2])-([5]). The algorithm solves the Vlasov equation 
in phase space and uses the so-called "splitting scheme" , 
a widely adopted strategy in numerical fluid d yna mics. 
Such a scheme, pioneered by Cheng and Knorr [20[, was 
first applied to the study of the Vlasov-Poisson equations 
in the electrostatic limit and then employed for a wide 
spectrum of problems [3] . For different values of the pair 
(Mo,e), which sets the widths of the initial water-bag 
profile, we performed a direct integration of the Vlasov 
system (j2j)-(j5j). After a transient, magnetization is shown 
to eventually attain a constant value, which corresponds 
to the QSS value observed in the HMF, discrete, frame- 
work. The asymptotic magnetizations are hence recorded 
when varying the initial condition. Results (stars) are re- 
ported in figure QJa) where Mqss is plotted as function 
of e. A comparison is drawn with the predictions of our 
theory (solid line) and with the outcome of N-body sim- 
ulation (squares) based on the Hamiltonian (pQ), finding 
an excellent agreement. This observation enables us to 
conclude that (i) the Vlasov equation governs the HMF 



dynamics for N — > oo both in the homogeneous and non 
homogeneous case; (ii) Lynden-Bell's violent relaxation 
theory allows for reliable predictions, including the tran- 
sition from magnetized to non-magnetized states. 

Deviations from the theory are detected near the tran- 
sition. This fact has a natural explaination and raises 
a number of fundamental questions related to the use 
of Vlasov simulations. As confirmed by the inspection of 
figure [2(b), close to the transition point, the entropy S of 
the Lynden-Bell coarse-grained distribution takes almost 
the same value when evaluated on the global maximum 
(solid line) or on the saddle point (dashed line). The en- 
tropy is hence substantially flat in this region, which in 
turn implies that there exists an extended basin of states 
accessible to the system. This interpretation is further 
validated by the inset of figure [2(a), where we show the 
probability distribution of Mqss computed via N-body 
simulation. The bell-shaped profile presents a clear peak, 
approximately close to the value predicted by our theory. 
Quite remarkably, the system can converge to final mag- 
netizations which are sensibly different from the expected 
value. Simulations based on the Vlasov code running at 
different resolutions (grid points) confirmed this scenario, 
highlighting a similar degree of variability. These findings 
point to the fact that in specific regions of the parame- 
ter space, Vlasov numerics needs to be carefully analyzed 
(see also Ref. [21]). Importantly, it is becoming nowadays 
crucial to step towards an "extended Vlasov theoretical 
model which enables to account for discreteness effects, 
by incorporating at least two particles correlations inter- 
action term. 
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FIG. 1: Panel (a): The magnetization in the QSS is plotted 
as function of energy, e, at Mo = 0.24. The solid line refers to 
the Lynden-Bell inspired theory. Stars (resp. squares) stand 
for Vlasov (resp. N-body) simulations. Inset: Probability 
distribution of Mqss computed via N-body simulation (the 
solid line is a Gaussian fit). Panel (b): Entropy S at the sta- 
tionary points, as function of energy, e: magnetized solution 
(solid line) and non-magnetized one (dashed line). 

Qualitatively, one can track the evolution of the sys- 
tem in phase space, both for the homogeneous and non 
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FIG. 2: Phase space snapshots for (M ,e) = (0.5,0.69). 



homogeneous cases. Results of the Vlasov integration 
are displayed in figure [2] for (Mo,e) = (0.5,0.69), where 
the system is shown to evolve towards a non magnetized 
QSS. The initial water-bag distribution splits into two 
large resonances, which persist asymptotically: the lat- 
ter acquire constant opposite velocities which are main- 
tained during the subsequent time evolution, in agree- 
ment with the findings of [ll| . The two bumps are there- 
fore an emergent property of the model, which is correctly 
reproduced by the Vlasov dynamics. For larger values 
of the initial magnetization (Mo > 0.89), while keeping 
e = 0.69, the system evolves towards an asymptotic mag- 
netized state, in agreement with the theory. In this case 
several resonances are rapidly developed which eventu- 
ally coaelesce giving rise to complex patterns in phase 
space. More quantitatively, one can compare the veloc- 
ity distributions resulting from, respectively, Vlasov and 
N-body simulations. The curves are diplayed in figure 
[3] (a),(b),(c) for various choices of the initial conditions 
in the magnetized region. The agreement is excellent, 
thus reinforcing our former conclusion about the validity 
of the Vlasov model. Finally, let us stress that, when 
e = 0.69, the two solutions (resp. magnetized and non 
magnetized) [ll| are associated to a practically indistin- 
guishible entropy level (see figure [3] (d)). As previously 
discussed, the system explores an almost flat entropy 
landscape and can be therefore be stuck in local traps, 
because of finite size effects. A pronounced variability of 
the measured Mqss is therefore to be expected. 

In this Letter, we have analyzed the emergence of QSS, 
a universal feature that occurs in systems with long-range 
interactions, for the specific case of the HMF model. By 
comparing numerical simulations and analytical predic- 
tions, we have been able to unambiguously demonstrate 
that the Vlasov model provides an accurate framework 
to address the study of the QSS. Working within the 
Vlasov context one can develop a fully predictive theo- 
retical approach, which is completely justified from first 
principles. Finally, and most important, results of con- 
ventional Vlasov codes are to be critically scrutinized, 
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FIG. 3: Symbols: velocity distributions computed via N- 
body simulations. Solid line: velocity distributions obtained 
through a direct integration of the Vlasov equation. Here 
e = 0.69 and M = 0.3 (a), M = 0.5 (b), M = 0.7 (c). 
Panel (d): Entropy at the stationary points as a function of 
the initial magnetization: the solid line refers to the global 
maximum, while the dotted line to the saddle point. 



especially in specific regions of parameters space close 
to transitions from homogeneous to non homogeneous 
states. 
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